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bJ) Abstract 

_ We review recent progress in numerical relativity simulations of black-hole 

^ ■ (BH) spacetimes. Following a brief summary of the methods employed in the 

0^ . modeling, we summarize the key results in three major areas of BH physics: 

QQ . (i) BHs as sources of gravitational waves (GWs), (ii) astrophysical systems 

CN i involving BHs, and (iii) BHs in high-energy physics. We conclude with a list 

C*^ ■ of the most urgent tasks for numerical relativity in these three areas. 
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1. Introduction 

For almost a century now after its discovery, the Schwarzschild solution 
describing a static, spherically symmetric vacuum spacetime in Einstein's 
theory of general relativity (GR) has been a valuable tool for a wealth of 
theoretical and experimental studies, including in particular weak-field tests 
of the theory. For more than half a century, the validity of the solution 
beyond the Schwarzschild radius r = 2GMj(? and the consequent existence 
of BHs was thought to be a mathematical artifact of the field equations. This 
picture has changed dramatically in more recent decades. Crucial events in 
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this sense were the mathematical discovery of the Kerr solution describing 
rotating BHs [l|, and the astrophysical X-ray and radio observations in the 
1960s and 1970s. Cosmic X-ray sources, first detected in 1963 |^, were 
soon associated with compact stellar objects and, in particular, with mass 
exchange from an ordinary star onto a compact object in binary systems 
[3|. This interpretation received further support when the X-ray source Cyg 
X-1 was identified with the supergiant BOIb star HD 226868 jj, |5|, which 
established Cyg X-1 as the first observed X-ray binary. Over the next few 
years, the explanation of these systems in terms of accretion onto compact 
stellar objects in binary systems became widely accepted (see |6[ for a review). 
Measurements of the binary period and radial velocity provide lower limits 
on the mass of the compact companion. The values obtained for some X-ray 
binaries comfortably exceed the upper mass limit of about 3 Mq for neutron 
stars. Recent investigations of Cyg X-1, for example, predict 8.7 ± 0.8 Mq 
for the compact member |7|. 

Also starting in the early 1960s, optical and radio observations identified 
"star like" objects with surprisingly large redshifts |8|, |9|. If interpreted as 
cosmological in origin, the redshifts implied enormous distances and, accord- 
ingly, luminosities. Over the next decade the cosmological origin of these 
"quasi-stellar" sources (or quasars) became clear, and accretion onto super- 
massive BHs (SMBHs) |l0Ull | was accepted as the most plausible explana- 



tion of their energetics [12|, |l3 



By now, observations of stellar dynamics 
near the centers of galaxies and iron Ka emission line profiles provide strong 
evidence that many galaxies harbor BHs at their centers IJ, ll5[. BH masses 
are generally closely correlated to the bulge velocity dispersion and luminos- 



ity |l6l. Il7l|. For the Milky Way, astrometric and radial velocity observations 



of short-period stars orbiting the galactic center [18|, |19| predict an unseen 
gravitational mass of 4.1 it 0.6 x 10^ Mq. 

Electromagnetic observations provide strong but indirect evidence for the 
existence of astrophysical BHs. A more direct way of probing the BH nature 
of astrophysical compact objects is offered by the newly emerging field of 
GW astronomy. First-generation ground-based laser-interferometric detec- 
tors (such as LIGO, VIRGO and GEO600) have reached design sensitivity, 
and advanced versions of these detectors will be operational within a few 
years |20|, |2l|, |22|. One of the most prominent sources of detectable GWs 



for these detectors is the inspiral and coalescence of stellar mass and, pos- 
sibly, intermediate-mass BH (IMBH) binaries. At lower frequencies, space- 



based laser interferometric observations [23| have the potential to comple- 



ment ground-based efforts with high signal-to-noise ratio observations of mas- 
sive BH binaries. Space-based detectors are guaranteed to observe galactic 
binaries containing white dwarfs and/or neutron stars, and they may also de- 
tect the inspiral of stellar compact objects into massive BHs. Due to the weak 
interaction of GWs with any type of matter, including the detectors, digging 
physical signals from the noisy data stream represents a daunting task and 
heavily relies on so-called matched filtering techniques 2J| . Matched filtering 
is a common choice to search for signals of known form in noisy data, and 
works by cross-correlating the actual "signal" (i.e., the detector's output) 
with a set of theoretical templates. Accurate modeling of the GW sources 
thus plays a vital role in maximizing the scientific output from these experi- 
mental efforts. 

Event rate estimates for IMBHs and stellar-mass BH binaries are very 
uncertain, but advanced Earth-based detectors are expected to observe sev- 
eral events per year (see 25|] for a review). Belczynski et al. noted that 
recent metallicity measurements require revisions in the population synthe- 
sis models used for event rate estimates, suggesting that stellar-mass BH-BH 
binaries may be the first GW sources to be detected 26|. The observations 
of two possible X-ray binary precursors of BH-BH binaries provide a much 
needed observational constraint on compact binaries containing at least one 
BH, suggesting again that BH-BH binary mergers may be more common 



than anticipated [27j. Population synthesis models have large uncertainties, 
related e.g. to the poorly known common envelope phase, and some models 
would necessarily be ruled out if Advanced LIGO does not detect BH binaries. 
In this sense, GW detectors are guaranteed to put constraints on compact 
binary formation in the near future. Schutz 28| also pointed out that coher- 
ent data analysis using three or more detectors may further increase event 
rates by factors of a few, making prospects of doing astrophysics with Earth- 
based detectors even brighter. Estimates for event rates of SMBH binaries 
detectable by space-based interferometers range from 0.1 to 1000s per year. 



29, 30, 31, 32 



with a most likely estimate of ~ 20 — 30 

Over the past decade, BH modeling has also attracted a great deal of 
interest in the context of high-energy physics. According to the gauge/gravity 
duality, gravitational physics in anti-de Sitter (AdS) backgrounds describes 



field theories on the boundary of that spacetime [33|, |3J]. According to this 
duality, BHs in AdS are dual to an equilibrium field theory, and therefore 
interacting BHs or BHs off-equilibrium may represent interesting physics 
from the dual point of view. In fact, there is some evidence that some of 



the physics behind quark-gluon plasma formation in heavy-ion colhsions is 
captured by this duahty, which on the gravitational side corresponds to shock 



wave collisions and BH formation in AdS |35 . 

Understanding BH spacetimes is also of interest for high-energy physics 
in connection with extra-dimensional scenarios. Solutions to the hierarchy 
problem in physics have been developed where all standard-model interac- 
tions are confined to a 3+1 dimensional brane embedded in a spacetime with 



large extra dimensions or an extra dimension with a warp factor |36l . l37j . In 
these models gravity propagates in the entire higher-dimensional spacetime 
(the "bulk"), and the fundamental Planck scale could be as low as 1 TeV. 
These energy scales are accessible to experiments such as the Large Hadron 
Collider (LHC). If the fundamental Planck scale is so low, Thome's hoop 
conjecture suggests the exciting possibility of BH production in high-energy 



parton-parton collisions [38|, |39|, i.e. a direct observable signature of the ex- 
istence of extra-dimensions. Accurate modeling of energy and momentum 
losses in the form of GWs during the collisions and determination of the 
BH formation cross section will be crucial for the analysis of data from such 
experiments (see e.g. (40|), and they should be amenable to classical calcu- 



lations in D-dimensional GR 41 



Finally, the understanding of BH dynamics under extreme conditions is 
fascinating from a mathematical-physics point of view. Questions such as 
BH stability, cosmic censorship, etcetera can only be addressed by solving 
the full nonlinear set of the field equations. 

For many of these scenarios, it is convenient to divide the coalescence of a 
BH binary in the framework of GR into three stages which mirror the different 
tools used for their modeling: (i) The inspiral phase, where the interaction 
between the individual holes is still sufficiently weak to justify the use of post- 
Newtonian (PN) techniques [42]; (ii) The final orbits, plunge and merger, 
which are governed by the strong-field regime of Einstein's equations and 
can be described only by fully numerical simulations; (iii) The ringdown of 



the remnant BH, which is amenable to a perturbative treatment [43|, |44, |45 

In this article we will focus on the second, strong-field stage, on the 
numerical tools dedicated to its study, and results obtained following the 



2005 breakthroughs [46|, |47|, |48| . As we shall see, however, a comprehensive 
understanding of BH dynamics requires a close interplay of numerical and 
analytical studies, and we will discuss in various cases the interface between 
numerical relativity and approximation techniques. 

Following a brief summary of the framework employed in numerical rel- 



ativity (Sec. |2]), in Sees. IMS] we will present the main results obtained from 
numerical studies of BHs in the context of GW detection, astrophysics and 
high-energy physics, respectively. We conclude in Sec. [6] with a discussion of 
future applications of numerical relativity. At the end of each section we will 
provide references to the literature for further reading. 

Notation: We shall be using geometrical units, such that the gravitational 
constant G = 1 and the speed of light c = 1. We denote spacetime indices 
. . . D — 1 by Greek letters, and spatial indices 1 ... D — 1 by Latin letters. 

2. Numerical modeling of black holes 

In Einstein's GR, spacetime is modelled as a D-dimensional manifold 
with a metric (7^^ of Lorentzian signature — h + + . . .. With the exception 
of Sec. O we will restrict our discussion to D = 4, i.e. three spatial and one 
time dimension. The metric is determined by the Einstein equations 

Gal3 = Rap — -^da/sR = TajS, (l) 

with Tai3 = in vacuum. While analytic solutions (as, for example, the 
Schwarzschild and Kerr solutions and the Friedmann-Lemaitre-Robertson- 
Walker solution) have been found for systems with high symmetry, solutions 
for dynamical configurations without special symmetries generally require 
the use of numerical tools. For this purpose the Einstein equations need 
to be cast as an initial value formulation such that the Cauchy-Kowalewski 
theorem guarantees a uniquely determined evolution of appropriately cho- 
sen initial data. In the case of Einstein's equations, it is not obvious that 
such a formulation can be obtained: a straightforward calculation shows that 
Eq. (]T]) contains second time derivatives of the metric components Qab, but 
not of Qta- This fact mirrors the coordinate or gauge freedom characteristic of 
GR, and implies the existence of so-called constraint equations which impose 
conditions on the metric components and their first derivatives on each spa- 
tial slice with constant coordinate time t. In the following we will discuss two 
approaches which provide an initial value formulation of the Einstein equa- 
tions and have led to successful, long-term stable numerical evolutions: the 
generalized harmonic (GH) formulation employed in Pretorius' breakthrough 
BH binary simulations J46| and the canonical Arnowitt-Deser-Misner (ADM) 



moving puncture breakthroughs 




split [49[, reformulated by York 50| , which forms the starting point for the 



ADM 3+1 Split: Spacetime is decomposed into a one-parameter family 
of three-dimensional spatial slices. We choose coordinates adapted to this 
formulation, such that the coordinate time t labels each slice and x* de- 
notes points on the slice. On each hypersurface, there exists a unique time- 
like, future pointing unit normal field n" which defines a projection operator 
±"^ = g'^p + n'^np onto the hypersurface. The geometry of the hypersur- 
face is completely determined by the three-metric or first fundamental form 
'jij = l.^iU^ jQ^u = Qij. The coordinate freedom is represented by the shift 
vector I3i = gu and the lapse function a = \/ gtt — P^Pm (with /9* = 7*^/3^), 
which relate spatial coordinates on (and measure separation in proper time 
between) different hypersurfaces. The projections 1.°" il.^ jG ap ■, -L'^iGa/j.n^, 
Gf^i^n^n" of the Einstein equations then lead to six evolution equations for 
the three- metric jij, three momentum constraints and the Hamiltonian con- 
straint: 

^tl^J = rdml^+l^id^r + lmAr-'iaK, (2) 

dtK,^ = rdmK,,+Krn,djr + K^Ar-DiD^a 

+a {TZij + KK,j - 2Ki^K"'j) , (3) 

n = n + K^-K^nK'^ = o, (4) 

M' = Dm{irm-Y"'K)=0, (5) 

where TZij, TZ and Di denote the Ricci tensor, Ricci scalar and covariant 
derivative associated with the three-metric, and the extrinsic curvature Kij 
has been introduced to obtain a first-order system in time. 

It turns out that this formulation of the Einstein equations is only weakly 



hyperbolic [51[, and therefore not suitable for long-term stable numerical evo- 
lutions. Motivated by these stability problems, several modifications of the 
ADM system have been investigated. The most popular is the Baumgarte- 



Shapiro-Shibata-Nakamura (BSSN) formulation [52, l53|, which rearranges 
the degrees of freedom via a split of the extrinsic curvature into trace and 
trace-free parts, a conformal transformation and promotion of the contracted 
conformal Christoffel symbols to the status of independent variables: 

K = Y'K,„ A, = e-'^ (Xi, - I^jK] , 

f * = 7"^"? t„„ = ^7'""7^'= {dmlnk + dnlkm - d^lmn) • (6) 



The resulting set of evolution and constraint equations is given in Eqs. (15)- 
(30) of SJ]. Subject to minor modifications, such as replacing the conformal 
factor by X = e"^'''' or adding the auxiliary constraint arising from the defini- 
tion of r* in Eq. ([6]), this formulation is employed in the present generation 



of moving puncture evolution codes (55|, |56|, l57|, l58|, l59l |60 



61 which fix the 



gauge via "1+log" slicing and the F-driver condition |5J, |62|, |63|, |64 . 

GH formulation: Harmonic coordinates are defined by the condition 

ga^dx^ = — Fq, = 0, where D = V^V^ represents the scalar wave opera- 
tor. These coordinates have played an important role in the analysis of the 



Cauchy problem in GR [65|, l66|, |67| . In harmonic coordinates, the Einstein 
equations take on a manifestly hyperbolic form, which allows for a general- 
ization to arbitrary coordinate or gauge choices 68|, |69|. The first step is to 



introduce four arbitrary source functions if" such that the coordinates obey 



F" = Dx' 



H", 



(7) 



and treat these functions as independent evolution variables. Then one con- 
siders the GH system 

-Rq/3 — V(qC/3) = 0, (8) 

where Ca = Ha + Ta- Eq. (^ is equivalent to the Einstein equations, subject 
to the validity of the constraint ([7]). In expanded form, the GH system is 
given by 

29"" g""^ {d^g^a.dxg,p - F„^«F^,a) , (9) 



g^^'d^d^g. 



-2V{aHp) 



and the constraints C^ = are preserved by virtue of the Bianchi identities 
provided that C^ and dtCa vanish on the initial hypersurface t = 0. A key 
ingredient in the numerical evolution of the GH system is the addition of a 
constraint damping term oc [25'^(Q,t/3) — gapt^\ C^ to the right hand side of 
Eq. ([9]) 70|. Here, t'^ is a non- vanishing timelike vector field. Finally, it 



is interesting to note that the Hamiltonian and momentum constraints are 



automatically satisfied on the initial hypersurface if Cq, = = dtCa [71 



The coordinates were determined in Pretorius' initial simulations by 



DH, 



<i 



a 



1 



a 



+ ^2n''d,Ht, Hi = 0, 



(10) 



and the evolution of the metric proceeds according to Eq. iQ with the con- 
straint damping term of Gundlach et al. l70| . For further discussion of gauge 
choices in the GHG system, see also 72|, |73|. 



The spectral code originally developed by the Caltech- Cornell group (see 



74,|75|] and references therein) employs a first-order version of the GH system 



7l[ with dual-coordinate frames and BH excision. Furthermore, the first- 
order GH formulation facilitates the specification of constraint-preserving 
boundary conditions 7l|. During the inspiral phase, they evolve the gauge 
functions H^ such that they remain constant in a frame comoving with the 
BHs, but switch to a dynamic evolution during the plunge and ringdown; see 
Sec. Ill in 76|. While more complex in structure, this framework enables 



their code to generate BH evolutions with exceptional accuracy |77l . 1 76 



Initial data: The construction of initial data requires solving the Hamilto- 
nian and momentum constraints. Here we only summarize the key concepts; 



we refer the reader to Cook's review article 781 for details. Most work on the 



construction of initial data is based on the York-Lichnerowicz split [79|, |50 



which involves a conformal transformation of metric and extrinsic curvature 
and separates the latter into trace and trace-free part. More recently, the 
thin-sandwich construction |80|, which replaces the extrinsic curvature by 
the time derivative of the metric, has become a popular alternative. In ei- 
ther case, the resulting elliptic equations simplify substantially under the 
assumption of conformal flatness and a spatially constant K. 

While this formalism provides a convenient method to solve the con- 
straint equations, we still need to ensure that the initial data represent a 
physically realistic system, typically two BHs with specific spins and mo- 
menta. This can be achieved by generalizing the Schwarzschild solution, 
which is obtained in the above framework in conformally fiat form with con- 
formal factor exp(0) = 1 + ^. A generalization to initial data of N BHs 

"'or 



81 



starting from rest is directly obtained by the construction of Misner 

Brill and Lindquist |82|. Remarkably, an analytic solution for the extrinsic 

curvature can still be obtained for BHs with initial linear momenta P„ and 



spins S„ 831]. By applying a compactification to the internal asymptotically 
fiat region, Brandt & Briigmann 84] arrived at the so-called puncture data 
construction, which is the method of choice for most numerical evolutions 
employing the BSSN formulation. In order to overcome an upper limit of 
~ 0.93 for the dimensionless spin of BHs in puncture type initial data 85| . 



Lovelace et al. (86| generated initial data based on a generalization of the 
Kerr-Schild solution 87|] and were thus able to evolve a BH binary with spin 
magnitude 0.95. 

In alternative to generalizing analytically known BH solutions, the pres- 



ence of horizons in the initial data can be encoded in the form of boundary 
conditions for the metric and extrinsic curvature 
the isolated horizon framework 
constructed in 



90 



891 as determined by 



88, 89 



Initial data along these lines have been 
and form the starting point of most of the numerical 
evolutions using the GH system (see e.g. 9l|, |75|). 

Diagnostics: Extracting physical information from numerical simulations 
of the Einstein equations is nontrivial for two reasons. First, the evolved 
variables are dependent on the coordinate conditions; second, it is often not 
possible to define local quantities familiar from Newtonian physics. The 
first difficulty requires the calculation of gauge-independent variables. The 
second difficulty is alleviated by the isolated horizon framework [90], which 



facilitates the calculation of BH mass and spin in the limit of isolated BHs; 
in practice, this framework is often applied when the BHs are farther apart 
and their interaction is considered negligible. Local properties of the BHs 



are encoded in their apparent horizon [93|, |9J, |95|, |96|. In particular, the 
irreducible mass can be expressed in terms of the apparent horizon area: 
Mirr = -^/Aah/IStt. In the limit of an isolated BH, the spin can be derived 
from the integration of the rotational Killing vectors over the horizon ac- 



cording to Eq. (8) of 961]. In practice it has been found that using flat-space 
rotational Killing vectors yields reasonable results 97|, but we also note the 
improved method to compute approximations to Killing vectors in 98| and 
an alternative method 99| based on Weyl scalars to compute quasi-local BH 
quantities without solving the Killing equation. It is important to bear in 
mind, however, the approximate nature of any spin calculation in BH binary 
simulations. 

In comparison, it is more straightforward to extract global quantities of 
the spacetime, most notably the total mass-energy Madm and the linear and 



angular momentum Pi, Ji |49|, |50[. These are given by surface integrals at 



92 



spatial infinity, see e.g. Eqs. (7.15), (7.56) and (7.63) in 

Arguably the most important diagnostic quantity in BH simulations is 
the GW signal. The most comrnon method to extract GWs is based on 
the Newman-Penrose formalism 100| and derives the complex Newman- 



Penrose scalar ' ^^ fr om contraction of the Weyl tensor with suitably chosen 



tetrad vectors 101 , Sec. HI A in 59 . It is often convenient to decom- 



pose \E'4 into multipoles ifjim using spherical harmonics of spin-weight —2: 
\E'4(t, r, 6^, 0) = ^;^2 -2^im(6', 0)'0«m(^5''^)- Contemporary numerical codes 
typically evaluate \l/4 at finite coordinate radius. This results in system- 



atic errors, due to ambiguities in the tetrad choice and the neglect of GW 
backscattering. While these errors can be reduced by extracting GWs at a 
sequence of rad ii and extrapolating to infinite radius (\E'4 asymptotically falls 
off ~ r~^ 102] ). a cleaner method is to calcul ate \ 1^4 at infin ity as provided 



by Cauchy- characteristic extraction methods |103l . Il04 Il05l ]. From \E'4, one 
straightforwardly obtains the energy, linear and angular momentum radiated 
in the form of GWs; see e.g. Eqs. (49)-(51) in 
is more common to work with the wave strain 



59| . In GW data analysis it 



h 



ihx 



dt' 



dt"^ 



4, 



(11) 



decomposed into multipoles: h{t,r,6,(j)) = ^^2~^^irn{G,4')him{t,r). In 
practice, after integrating \E'4 twice in time the signal can be severely af- 
fected by nonlinear drifts, which can be controlled t o a s i gnifi cant extent by 
performing the integration in the frequency domain 106l Il07 . 

Alternatively, GWs can be extracted by viewing the metric in the far field 
regime as a perturbation of th e Sc hwarzschild spacetime and employing the 



formalism of Regge, Wheeler |108| and Zerilli |l09l . Ill0| . One thus obtains 



the Regge- Wheeler- Moncrief and Zerilli-Moncrief master functions Q^^, Qtrai 
which can be converted i nto the multipolar components of the GW strain h 



according to Eq. (49) of [111 



For a comprehensive summary of the numerical framework for evolving 
Einstein's field equations, we refer the reader to the books by Alcubierre 51 



and Baumgarte and Shapiro |112J . Further details, including mathematical 



aspects of the equa tion s and cha racteristic techniques, can be found in the 



review articles 113, 92, 114, 115 



3. Gravitational wave physics 

BH binary systems represent one of the most promising sources of de- 
tectable GWs. The parameters of a BH binary are commonly divided into in- 



trinsic and extrinsic |116| . Intrinsic parameters characterize physical proper- 
ties of the system, such as the total mass M, the mass ratio q = M2/M1 < 1, 
the individual BH spins Si, S2 and the orbital eccentricity. Extrinsic param- 
eters such as sky position, distance, orbital inclination angle, arrival time and 
initial phase of the wave, in contrast, depend on the source location relative 
to the observer, and do not directly enter the GW source modeling process. 



10 



The majority of numerical studies of BH binary spacetimes performed 
to date have focussed on comparable mass-ratios q > 1/10 and moderate 
spin magnitudes Xi — |Sj|/Mf < 0.8, with particular emphasis on spins 
(anti-)aligned with the orbital angular momentum. We note, however, the 
following explorations outside this "best charted" subset of the parameter 
sp ace. Circu lar binaries with mass ratios up to g = 1/100 have been studied 
comparisons for this value of q with fully perturba tive calcu- 

who 



117, 118 



119 



m 

lations have been obtained in the limit of head-on collisions in 
report a discrepancy of ~ 7 % in GW energy and momentum, probably in 
large part due to the discretization error of the numerical simulations. The 
first direct comparison of numerical results with leading-order self-force pre- 
diction was studied in the context of periastron advance in BH b inar ies in 
BH binaries with spins 



120 



0.9 have been evolved in 121, 85, 122 



124, 61 



but exceeding an apparent barrier ofy ~ 0.93jl23| requires departure from 

The present maximum is 



the conformal flatness assumption 

Xi,2 = 0.95 for aligned spins evolved for 12.5 orbits by Lovelace et al. 

Fina lly, w e note that emission of GWs efficiently circularizes BH binary or- 



86 



bits |125| . BH binaries are therefore expected to have vanishing eccentricity 
by the time they enter the frequency window of ground-based detectors, 
and for this reason most work on GW source modeling has focused on the 
quasi-circular limit. The derivation of BH momenta, including a small ra- 
dial component, for quasi-circular initial data in numerical relativity is based 
either on integrating the PN equations of rnotion 1261. 127l| or iter ative proce- 
dures using several numerical simulations 128l . l75l . Il29l . Il30l . Il3l| . The decay 
of eccentricity during the inspiral phase as well as periastron advance has 

For numerical studies of 



129 



been measured in numerical simulations in 
BH binaries with significant eccentricity as well as arguments why eccentric 
binaries may represent relevant sources of GWs after all, especially in the 
more extreme mass ratio reg i me a nd for space-based detectors, we refer the 



reader to 132, 133, 134, 135, 136 



Gravitational waveforms from BH binaries: For illustration, we dis- 
play in Fig. [1] the real parts of the /i22 and /las multipoles of the GW strain 
h obtained for the inspir al an d coalescence of a binary of nonspinning BHs 



with mass ratio q = 1/4 [l37|. In the course of the inspiral, both amplitude 
and frequency of the GW signal gradually increase. Close to merger, defined 
as formation of a common apparent horizon, the GW amplitude reaches a 
maximum. Eventually it drops exponentially as the merged hole rings down 
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Figure 1: Real part of the GW multipoles /122 and ^133 obtained from the inspiral and 
merger of a nonspinning binary with mass-ratio q = 1/4. The dashed red curve represents 
PN predictions matched to the numerical (solid black) signal. The vertical dotted lines 
mark the time when the frequency of the (Z, m) = (2, 2) multipole has reached Mo; = 0.1. 



to a stationary Kerr state: see e.g. Fig. 18 in 91 . 

The qualitative features of the GW signal emitted by different types of 
binaries can be summarized as follows, (i) The inspiral of two nonspin- 
ning, equal-mass BHs is the case most intensively studied by numerical 
relativity. The wave signal is dominated by the (/, m) = (2, 2) multipo- 
lar component, which carries > 98 % of the total radiated energy. The 
longest numerical simulation by Scheel et al. [1381 ] covers 16 orbits and re- 
ports a ratio of final to initial BH mass Mf/M = 0.95162 ± 0.00002 and 
a final spin Sf/Mj = 0.68646 ± 0.00004. (ii ) As the mass-ratio decreases, 
higher-order modes become more prominent |l39l | and the fract ional energy 
in (/,m) = (2,2) drops to about 68 % in the limit q ^ [m^. Motivated 
by the stron g rel ation of phase and frequency of different (/,m) multipoles. 
Baker et al. |l41| view the (/, m) radiation modes as being generated by the 



corresponding momenta of an implicit rotating source; see also |142| . (iii) 



GW multipoles with odd I are suppressed in the equal-mass limit due to 
symmetry, (iv) BH binaries with spins aligned with the orbital angular mo- 
mentum emit stronger GW signals due to an effect sometimes referred to as 
orbital hang-up: equal- mass binaries with Xi = X2 = 0.757 radiate about 
twice as much energy and angular momentum when compared to the non- 



spinning case [143(1 . For anti-aligned spins the GW energy decreases by a 



12 



similar factor, (v) Orbital eccentric ity i ntroduces a nonmonotonic behavior 



of the GW frequency; cf. Fig. 6 in 133[. (vi) Spin-spin and spin-orbit cou 



plings cause a precession of the orbital plane in the case of BH binaries with 
spins that are not aligned or anti-aligned with the orbital angular momen- 
tum. This precession manifests itself in a modulation of the GW amplitude 
emitted in a fixed angular direction: cf. Fig. 1 in 144 . 



Tests of general relativity: Numerical simulations of binary BH mergers 
provide an opportunity to study the nonlinearities of the theory and to test 
the Kerr nature of astrophysical BHs. If GR is the correct theory of gravi- 
tation, all BHs in the Universe should be described by the Kerr solution [l|, 
which depends on two parameters: the mass M and the dimensionless spin 
X- The remnant of a binary BH merger settles down to the Kerr solution 
by emitting gravitational radiation at characteristic (complex) quasinormal 
frequencies Munim that depend only on x- Here (l,m) are the usual angular 
indices, and n is the "overtone number": modes with small n have a longer 



damping time and should dominate the radiation |145l . 145 



The measurement of the real and imaginary part of a single quasinormal 
mode contains, in principle, enough information to determine the mass and 
spin of the BH. Together with the accurate determination of the individual 
BH masses from the inspiral phase this already allows for tests of GR, by 
testing the GR prediction for mass loss during inspiral and merger and at 
the same time the strong field dynamics of GR. The measurement of any 
other frequency or damping tim e the n provides a further interesting test of 
the Kerr nature of the final BH 146| . The feasibility of these tests depends 



on the characteristi cs o f a gi ven GW detector and on the relative excitation 
of the modes (see |147l . Il48j for detailed studies). Quantifying the excita- 
tion of quasi norm al modes for generic initial data is a formidable technical 
problem 4^ Il49| , but numerical merger simulation s all ow sensible estimates 



of the relative quasinormal mode amplitudes |9l|, Il39| . The relative mode 
amplitude depends on the binary parameters, and therefore, in principle, the 
measurement of a multi-mode signal with single or (better) multiple GW 
detectors can be used to measure the binary mass ratio or the inclination of 
the final BH spin with respect to the line of sight. 

We also remark that the increased signal-to-noise ratio coming from the 
merger has been shown to improve the bounds on alternative theories of 
gravity that could come from observations of the inspiral only. Keppel and 
Ajith, for example, discussed this possibility in the context of massive gravi- 
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ton theo ries, where the graviton mass would modify the dispersion relation 



of GWs [150|. Formulations of the evolu tion equations in alternative theories 



of gravity are in their infancy |l5ll |152| and the numerical e xplo ration of BH 



binaries in alternative theories of gravity has just started [153[. This is an 
interesting line of research, and in the near future we may have more concrete 
ways to quantify deviations from GR in strong-field mergers. 

GW template banks: The key challenge in GW data analysis is to ac- 
curately predict and isolate the features of gravitational waveforms in the 
data stream from GW detectors. This is necessary in order to (i) detect 
the presence of a GW signal of astrophysical origin, and (ii) determine the 
parameters of the emitting source. These goals can be achieved by matched 
filtering, i.e. by cross-correlating the data stream s (composed of detector 
noise n plus a potential GW signal h) with a bank of theoretical waveform 
templates hx., where the indices Aj [i = 1, . . . ,P) denote the P intrinsic and 
extrinsic source parameters. For this purpose it is imperative to construct a 
large set of accurate numerical waveforms which cover the relevant parameter 
space. This task far exceeds the capacity of purely numerical methods. If we 
assume a seven- dimensional parameter space - one mass ratio and three spin 
parameters per hole - a mere ten waveforms per parameter would correspond 
to 10^ waveforms. Furthermore, using binary waveform models including the 
long inspiral phase plus merger optimizes the signal-to-noise ratio in GW 
observations, and the inclusion of merger and ringdown has been shown to 
provi de s ignificant improvements in source localization and distance calcula- 
tion 154l ]. Because current numerical simulations cover only about 30 cycles 
of the inspiral, the generation of complete waveforms requires the combined 
use of PN and numerical methods. The construction of such GW template 
banks currently proceeds along either of the following two paths. 



(i) The effective- one- body (EOB) approach |l55l . Il56l | maps the dynamics 
of the two-body problem in GR into the motion of a particle in an ef- 
fective metric. The components of the effective metric are currently de- 
termined to 3PN order. The EOB method improves upon this model by 
using additional pseudo-PN terms of higher ordeiij which are not derived 
from PN expressions, but calibrated via comparison with numerical results 



^Comparing an EOB model without such additional calibration against numerical rel- 
ativity waveforms has been found to result in an accumulated phase di fferen ce at merger 
of 3.6 rad for nonspinning binaries with mass ratios q between 1 and 6 |157| . 
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1571 . 11581 . 11591 . 11601 . 11611 . 1162 



Further improvements come from using a re- 
summed version of the PN expanded res ults a nd from modehng nonadiabatic 
effects in the inspirah see e.g. Sec. IV in 163| . The inspiral-plunge waveform 



resulting from this construction is then matched to a merger-ringdown sig- 
nal composed of a superposition of quasinormal oscillation modes of a Kerr 
BH. The total number of free parameters varies between the individual EOB 
models currently investigated, but in all cases, the general strategy is to cal- 
ibrate these parameters by comparison with a finite number of numerical 
simulations. 

(ii) Phenomenological waveform templates are based on hybrid waveforms: 
the early inspiral is modeled by PN techniques, and the resulting signal is 
matched (within a specified window) to a numerical waveform describing the 
last orbits, merger and ringdown. The resulting set of hybrid waveforms 
are then approximated by a model containing a number of phenomenological 
parameters which are mapped to the physical parameters of the binaries, 
such as the mass ratio q. T he p h enomenol ogical models initially developed 

and then extended to the case 
are based on 



for nonspinning binaries in |164| . Il65l |166 



of spins aligned with the orbital angular momentum |167 



hybrid waveforms matched i n th e time domain, whereas the more recent 



study by Santamaria et al. [l68f | performs the matching in the frequency 



domain. In spite of such differences in their construction, the final result of 
all these phenomenological models are closed-form analytic expressions for 
the waveforms in the frequency domain. An exploration of phenomenological 
models for equal-mass binaries with spin precession is given in 



169 



Both types of template banks are employed in the analysis of GW data. 
The recent search for GWs from binary BH inspiral, merger and ringdown 
170| used phenomenological models for injection, and an EOB model for 
injection and for the search templates. This search constrained the rate of 
mergers of binaries with individual masses in the range [19, 28] Mq to be no 
more than 2.0 Mpc~'^Myr~^ at 90% confidence. Numerical w avefo rms have 



also been used inside the GW community-wide Ninja project |17l| to study 
the sensitivity of existing GW search algorithms used in the analysis of obser- 
vational data. For this purpose, numerical relativity waveforms were injected 
into a simulated data stream designed to mimic ground-based detector noise. 
The algorithms detected a significant fraction of the injections, but likely re- 
quire furthe r deyelopment, in particular for the purpose of measuring source 



paramet ers jl72l . Il73l . Il74l |. A further community wide effort, the NRAR 



project [175|, is dedicated to a systematic exploration of the complete BH 
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binary parameter space to develop optimally-calibrated template families to 
be used in GW data analysis. 

One of the most exciting prospects of GW detection is the idea of looking 
for electromagnetic counterparts to the GW signal. While a network oi Earth- 
base d de tectors is required to localize a GW source by triangulation (see 
e.g. |l76| ). a single space-based detectors like LISA may suffice to localize 
low-redshift sources and determine their luminosity distance. The inclusion 
of merger and/or of higher multipolar components of the radiation in the 
waveform models has been shown to provide s ignificant improvements in 



177. 154 . For this reason 



source localization and distance determination 
it is very important to construct phenomenological models including higher 
multipoles of the full inspiral /merger /ringdo wn w aveform. The EOB model 
has recently been extended in this direction 157 . 

Accuracy requirements: Clearly, the two approaches for template con- 
struction mentioned above require accurate numerical waveforms. Quantify- 
ing these accuracy requirements is a nontrivial task. Numerical uncertainties 
due to finite resolution are directly tested by convergence analysis, and the 
error incurred by extracting GWs at finite radius can be estimated using 



extrapolation of results from various radii to infinity |178| . The uncertainties 



in phase and amplitude thus derived, however, dep end on the alignment of 
the waveforms in phase and time; see e.g. Fig. 2 in 137|. Because offsets in 
phase and time are free parameters in GW data analysis, it is more conve- 
nient to measure accuracy and agreement of different waveforms in terms of 
quantities which take such alignment into account by construction. Such a 
measure is obtained from the inner product betw een two w aveforms h{t) and 



g{t) used in the theory of parameter estimation 179l Il76 



{h,g) = me 



SnU) 



df, 



(12) 



where the tilde and asterisk denote Fourier transform and complex conju- 
gate, respectively, and SV(/) is the one-sided power spectral density of the 
detector strain noise 180|. In practice, this inner product is to be under- 
stood as maximized over constant offsets in time and phase (Ato and A0o) 
between the two waveformqj. The signal-to-noise ratio that would be ob- 
tained for a physical signal he and a model waveform hm is then given by 



^For detection of GW events, one often considers the effectualness of waveform models 
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{he\hm) /\\hm\\- It IS related to the opti mal s ignal-to-noise ratio p for a 



Pr, 

perfect model waveform by the mismatch M. |176 



p^ = (1 - M)p = (1 - M){K\h,)/\\K 



(13) 



To leading order in ||5/i||, this definition of the mismatch implies 18ll . Il82 
^ ^ {5h\5h) - {5h\\\5h\\) ^ {5h\5h) ^^^^ 



2{he\K) 



2{K\he)' 



where 6h\\ = {6h\he) / {he\he) is that part of the waveform error parallel to he- 
Based on these definitions, accuracy requirements are obtained as fol- 
lows. Two wav eforms di ffering by 6h will be indistinguishable for a detector 



if {6h\6h) < 1 |l8lL Il63| . such that acceptable waveform models to be used 
in parameter estimation must lie within a sphere of unit radius of the exact 
waveform. For the purpose of event detection, we note that the fractional 
loss in detection of GW signals due to imperfect templates is proportional to 
the third power of the reduction in the SNR, i.e. ~ 3A^. A value of 10% is 
typically deemed acceptable, corresponding to a mismatch Ai = 3% . . . 3.5%. 
For detection efficiency, there arises, however, a complication due to the dis- 
crete nature of the template bank. As discussed in detail in Sec. II B of 



181 



the effective mismatch in a real GW search is given by a contribution 
due to the actual modeling procedure, plus a term which accounts for the 
discrete spacing in the parameter space: see their Fig. 3 and Eq. (17). Tem- 



plate banks adopted in LIGO sear ches use 0.3 for the latter 170l . Il83| and. 



Lindblom at al. |181| recommend a maximum mismatch 
In term s of the right-hand side of Eq. flT^ , dubbed inac- 



m consequence, 

A^max = 0.005. 

curacy function in Ref. [163[, we summarize the accuracy requirements as 



(cf. Eq. (10) in Ref. |l82| ) 



\m 

WhW 



< 



1/p 



V2Mr 



for parameter estimation, 
for detection. 



(15) 



Accuracy measurements employing the mismatch or \\Sh\\ have been used 
to study numerical, semi-analytic and hybrid waveforms. According to the 



(as opposed to faithfulness) and in those cases also maximizes the inner product over the 
source parameters. 
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Samurai project 184l |. relative mismatches of the I = m = 2 mode be- 
tween a variety of purely numerical waveforms for bina ries with total masqj 
M > 60 Mq are better than 10~^. Hannam et al. 185| found errors in 
the construction of hybrid waveforms to be dominated by PN contributions. 
Numerical waveforms should be long enough to enable matching at lower 
frequencies, where PN approximations are more accurate: ~ 3 orbits before 
merger for the equal-mass, nonspinning case , and ~ 10 orbits for binaries 
with spins x = 0-5- MacDonald et al. 182| report significantly more de- 
manding requirements of ~ 30 orbits for nonspinning, equal-mass binaries. 



Similarly, Boyle |186| concludes that longer NR waveforms or more accurate 
PN predictions will be required. The discrepancy largely arises from the use 
of more stringent acc uracy thresholds: A^max = 0.005 in 182| (instead of 
-A^max = 0.03 in 185|) for de tecti on, and p = 40 for parameter estimation. 



Most recently, Ohme et al. |187l ] have argued that maximising the match 
over intrinsic source parameters as well as phase and time shifts may reduce 
length requirements on numerical waveforms. 

A comparison of phenomenological and EOB template banks bas ed on 
the inaccuracy function has been performed by Damour et al. 163|. For 



advanced ground-based detectors, they conclude that current phenomeno- 
logical models deviate from the EOB (considered as a target model) beyond 
acceptable thresholds over a wide range of the parameter space, for both, 
detection and parameter estimation. 

Furt her reading on BH simulations in the context of GW physics is avail- 



able in 188, 189, 190, 191 



4. Numerical relativity and astrophysics 

Astrophysical observations have provided the strongest evidence to date 
that BHs exist and play an important role in many physical processes in our 
Universe. Astrophysical observations of X-ray binaries or quasars, as men- 
tioned in Sec. [H have a good deal to tell us about BHs. In recent years, 
numerical relativity simulations of BHs have also deepened our insight into 
a variety of astrophysical systems. We will first summarize what we have 
learned from "traditional" electromagnetic observations about the popula- 



■^For significantly lower masses, the last 10 orbits and merger signal lie outside the 
maximum sensitivity range of the eurrent generation of ground-based GW deteetors. 



tions of BHs we can expect to exist in the Universe. Then we will discuss 
some astrophysical implications of numerical relativity. 

Expected BH populations: We have already mentioned the two main 
classes of BHs identified by astrophysical observations, (i) So lar-m ass BHs 
with M ~ 5 — 20 Mq are usually found in X-ray binaries |l92| and are 
formed as the end-product of the evolution of massive stars, (ii) SMBHs 
with M ~ 10^ — 10^'^ Mq are believed to reside in most Active Galactic 
Nuclei (AGN) [13| . The assembly of these SMBHs is likely to result from a 
combination of BH mergers and accretion of surrounding matter over cosmo- 
logical timescales. Evidence for a third class of BHs with M ~ 10^ — 10^ Mq 
(IMBHs) is tentative at present flOJ. Eol. E95|. 

GW frequencies are inversely proportional to the system's total mass, and 
therefore BH masses play a crucial role in binary BH detectability. Earth- 
based detectors, such as LIGO and Virgo, have an optimal sensitivity band 
corresponding to stellar-mass BHs and IMBHs, while the planned space- 
based detector LISA is most sensitive to high-mass IMBHs and SMBHs. As 
discussed in the previous section, the expected GW pattern from BH binaries 
depends on the masses and spins of the binary members. We shall further see 
below that the gravitational recoil imparted upon the remnants of coalescing 
binaries strongly depends on spins and mass ratio. So, what do "traditional" 
electromagnetic observations tell us about BH masses and spins? 

Current spin estimates for accreting, stellar mass BHs are usually ob- 
tained by modeling the shape of their X-ray spectrum or by analyzing the 
skew in Fe Ka emission lines. The resulting estimates are all, to some extent, 
model-dependent, but frequently yield moderate values x ~ 0- 1 • • • 0.8 and, 
for some cases, values close to the Kerr limit x = Ij see e.g. 45|, Il96| and ref- 
erences therein. Theoretical arguments and observations suggest that stellar- 
mass BHs in binaries retain the spin they had at birth: neither accretion nor 
angu lar mome ntum extraction are likely to change significantly their mass or 



spin |l97l . Il98| . For SMBHs, spin estimates relying on the sh ape of the Fe Ka 



line yield values ranging from moderat e y ~ 0.6 . . . 0.8 (e.g. 1991 . l200l 1201 



to 



near-critical spins (e.g. 202l . |203| . l204j . |205|). SMBHs are expected to grow 



by a combination of mergers and accretion. Their s pin depends se nsitively 
on the details of these processes and of their growth 206l 12071 . |208| . SMBH 
assembly via BH mergers does not appear to be able to account for the large 
observed spin values; this may indicate that accretion dominates over merg- 
In any case, BH spins encode the history of their formation, and 



ers 



209 
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it would be extremely useful to have detailed knowledge of the BH spin dis- 
tribution function. Ultimately, such a detailed c ensus of BH parameters is 
one of the key targets of future GW observations 210l . l21ll ]. For the purpose 



of numerical modeling of BHs, present measurements indicate that all spin 
magnitudes (0 < x < 1) should be covered. 

We have already noted that stellar-mass BHs typically have masses in 
the range [5,20] Mq. Unfortunately, there are no confirmed observations 
of stellar-mass binary BH candidates, so estimates of mass ratios must rely 
on theoretical models. For stellar-mass binaries, population synthesis codes 



suggest that q should always be quite close to unity |198| . Measurements of 



SMBH masses are obtained by observations of stellar motion near galactic 
centers, as mentioned in Sec. [H as well as motion of gas discs J212l| . applica- 
tions of the virial theorem to the velocity dispers ion o f stars IJ] and rever- 



beration mappings applied to more distant AGNs [2 131 ]. An exhaustive list of 
galaxies with SMBH mass measurem ents in the range M ~ 10^ — 10^ Mq is 



presented by Graham et al. |214l . |215| ] . The impact of different SMBH assem- 
bly models on the mass and mass ratio distribution of detectable binaries has 
been discussed by various authors. The general consensus is that m ass ratios 



g < 1/10 (and down to g ~ 10 ^) should be common 3l|, |216| . 12171 ]. In con 



trast with the case of stellar -mass BHs, there i s by now some observational 
evidence for SMBH binaries [iisl, [219I, O, H, [iij • These observations are 
not sufficient to tightly constrain SMBH merger rates, but they are at least 
broadly consistent with merger-tree models predicting tens to hund reds of 
events during the typical lifetime of a space-based GW detector 223 . 



A more speculative kind of source for Earth-based detectors consists of 
the intermediate mass ratio inspiral of a compact object (neutron star or 
BH) into an IMBH. Another promising source for advanced Earth-based 
interferometers (albeit with highly uncertain event rates) are IMBH-IMBH 
mergers. The initial inspiral of these binaries could be detected via space- 
based interferometers, while the ringdown phase is in the optimal bandwidth 
for second- and third-generation detectors such as Advanced LIGO and ET, 



that could therefore be used for "follow-up" ringdown searches [2241 . |225 



As in the case of BH spins, the observations imply that numerical rela- 
tivity needs to cover a wide range in q. This can be done most efficiently by 
bridging the gap between numerical studies and t he perturbative modeling 
of extreme-mass-ratio binaries (e.g. 2261 . 12271 1228| ). 

From this discussion, it is clear that a detailed knowledge of the spin evo- 
lution as well as the generation of gravitational recoil in BH binary mergers 
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is important for understanding the cosmological evolution of SMBHs over 
cosmological times. We will discuss these effects in turn. 

Black-hole spins resulting from mergers: Prior to the 2005 numerical 
relativity breakthroughs, it was known that "minor mergers" {q < 10) of 
a large, rotating BH with an i sotropic distribution of small objects would 
tend to spin down the hole 229l |. Numerical merger simulations showed that 
the merger of comparable-mass, nonspinning BHs leads to a final Kerr BH 
with spin parameter a/M = 0.69. The simulations were followed by the 
development of several models to predict the spin of merger remnants as a 
function of the binary parameters for generic mass ratios and spins. 

The first studies focussed on performing numerical evolutions of the last 
few orbits of BH binaries, and used the results to calibrate formulae that 
map the initial binary pararn eters to values for the spin of the merged hole 



97, 139, 121, 230, 231, 232 



It became clear, however, that the binary 
inspiral up to the last orbits has the potential to significantly affect the spin 
distribution (e.g. 233[) and therefore should be included, for example via PN 



modeling, in the de rivat i on of maps fr o m in itial parameters to the merger 



remnant properties [2341 . |235| . |236| . 12371 . |238| . Predictions for the final spin 



based on the extrapola tion of te st-particle calculations to finite mass ratios 
have been developed in 158l . l239t |. and show remarkably good agreement with 
numerical results in the comparable-mass regime. A comprehensive review 
of all spin formul ae is beyond the scope of this article, but w e refer the reader 
to the review in 240l | and the discussion in Sec. V of 238| . If we assume an 
ensemble of BH binaries with initially r andomly orien ted spins, the final s pin 
distribution is peaked around x/ ~ 0.7 J208l . l234l . l237l |. Campanelh et al. [9 71 . 
24 1| reported spin flips by 34 — 103° with respect to the initial spin direction 
of the larger hole; spin flips of th is magnitude coul d pro vide an explanation 



for X-shaped radio sources [242 



238 demonstrated that 



Kesden et al. 

spin precession over the course of a long inspiral of thousands of orbits tends 
to align (antialign) the binary BH spins with each other if the spin of the 
more massive BH is initially partially aligned (antialigned) with the orbital 
angular momentum, thus increasing (decreasing) the average final spin. Spin 
precession is stronger for comparable masses, and it could produce significant 
spin alignment before merger for both SMBHs and stellar-mass BH binaries. 

Gravitational recoil: One of the most spectacular results obtained from 
numerical BH binary simulations is the quantitative prediction of the mag- 
nitude of gravitational recoils (or kicks). In the 1960s it was realized that 
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the emission of linear momentu m via GW s must impart a kick on the source 



due to momentum conservation [2431 . l244j . In the 1980s, the kick generated 



in compact binary inspirals and plu nges was s tudi ed in the framework of PN 



theory and BH perturbation theory |245l . l246l . l247l |. However, the astrophysi- 



cal relevance of the gravitational recoil following BH binary mergers remained 
an open question until the recent numerical relativity breakthroughs. 

For the case of an equal-mass, nonspinning BH binary, the net linear 
momentum emitted in GWs vanishes due to symmetry. Nonzero recoils are 
therefore only generated in systems where this symmetry is broken through 
(i) a mass ratio g < 1 (equivalently, a symmetric mass-ratio parameter rj = 
q/{l + qY < 1/4) or (ii) nonvanishing spins. 

For zero spins, early numerical studies in the range q = [1,1/4] found 



that the kick velocity is well approximated by J57l . l248l . |249| 



■^-'kick 



1.2 X 10VVl-4r/(l - 0.93r/) km/s. 



(16) 



This result (whose functional form is inspired by Fitchett's analytical wor k 
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was confirmed by simu lations for smaller mass ratios (g = 1/10, [140|) 
and by analytical work [250|. The maximal recoil obtained from Eq. ( IT6l) is 
■^max = 175.2 ±11 km/s for q = 0.36 ±0.03. Sm all or bital eccentricity e < 0.1 



should introduce corrections proportional to e [251 



Spinning binaries are characterized by seven free parameters (the mass 
ratio plus three components for each BH spin), and numerical studies in- 
evitably focussed on subsets of the parameter space. It soon became clear 
that the spin interaction dominates over mass ratio effects. The first studies 
for equal-mass binaries with spins parallel to the orbital angular momentum 
revealed kicks of several hundreds km/s, wit h an extra polated maximum of 
~ 500 km/s for extremal spin magnitudes 252l |253J . Shortly thereafter, 
the discovery of the so-called superkicks marked one of the most surprising 
outcomes of numerical relativity: binaries for which the spins are perpen- 
dicular to the orbital angular momentum and anti-aligned with each other 
can generate kicks of thousands of km/s, with an ex trapolated ma ximum 



4000 km/s for near-extremal spin magnitudes 24ll . l254l . |255[ |. Most 



recently, BH binaries with spins partially aligned with the orbital angular 
momentum but whose spin projections into the orbital plane still correspond 
to the superkick configurations have been foun d to result in even larger max- 
imum recoils of up to ~ 4 900 km/s J256l . l257l | , so-called "hang-up kicks" . 

Many galaxies harbor SMBHs at their centers, and galaxy mergers should 
generally lead to the merger of their central BHs 258| . Typical escape veloc- 
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ities range from ~10 km/s for dwarf galaxies up to ~ 1000 km/s for giant 



elliptic galaxies |259| . Large kicks would displace or eject the merged hole 
from its host, with possibly observable consequences: a softening of the stellar 
density gradient in the galactic nucleus, off-center radio-loud active galactic 
nuclei, off-nuclear X-ray sources in nearby galaxies and the generation of elec- 
trornagnetic s i gnals via inte r action of the BH with its gaseous environment 



2601 . I26ll . I262L I259L l263l . 1264 l265l . I266L l267l . l268l . l269j . BH ejection represents 



a potential obstacle for BH growth via merger, and thus puts constraints on 
merger-history mod els, which mus t be able to explain the assembly of SMBHs 



by redshifts z > 6 |26ll . l270l |271| . Observed redshifts of broad-line relative 



to narrow-line regions in quasar spectra may be interprete d as a smoking 



gun of BH ejection due to gravitational recoil |272l . |222| . |273| . but alternative 



interpretations (such as a BH binary, or s uperposed era ission regions from 
two interacting galaxies) are possible 274 |275| . l276l . |277 . 



The interest in astrophysical conse quences o f large reco ils le d to nume rical 



studies of the BH parameter space (278i |27a |280|, l28l|, Ififl |282|, l85L [283, 



284 |285| . Phenomenological kick formulas inspired by PN studies |144| | and 
similar to Eq. f lT^ were proposed to map the input parameters of a given 
binary configuration to the final kick velocity. Available numerical results 
span mass ratios in the range q > 1/ 10 a nd spin m agnitudes \x\ < 0.9, and 
are well described by Eq. (2) of Ref. |24l| (see also O, [286J). 

The apparent incompatibility between large recoils and the existence of 
BHs at galactic centers can be resolved by mechanisms that would align the 
individual BH spins with the orbital angular momentum of the binary. One 
such mechanism are gas torques in the so-called "wet" (gas-rich) me rgers , 
that would produce partial alignment "early on" in the inspiral phase 287 . 
If partial alignment in gas-rich mergers is the norm, as suggested also by 
the spin measurements discussed above, PN spin effects will lead to further 
alignment of the spins with the orbital angu lar mome ntum, significantly re- 



ducing the typical values of recoil velocities [233l . l288l |. This mechanism has 



even been shown capable of efficient suppressio n of large recoil velocities in 
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the case of the above mentioned hang-up kicks 

Merger simulations with matter and the Blandford- Znajek effect: 

Many numerical relativity groups are presently working on binary BH simula- 
tions in the presence of matter. Most of this work is trying to understand the 
signature of electromagnetic counterparts to binary BH m erger s: for exam- 



ple, such a signature could be produced by merger events |290l | or recoiling 
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BHs 291l | "shocking" the surrounding accretion disks. The Georgia Tech 



group presented the first simulations in full GR of equal-mass, spinning BHs 
merging in a gas cloud. They found that shocks, accretion and relativistic 
beaming can produce electromagnetic signatures co rrela ted with GWs, espe- 
cially when spins are aligned with the orbital axis 292|. Later work by the 



group focussed on hot, radiatively inefficient accretion flows. In this case, 
BH binaries exhibit a flare followed by a sudden drop in luminosity associ- 
ated with the m erger, an d quasi-periodic oscillations correlated with the GWs 



during inspiral [293|, l294j . The Urbana group simulated equal-mass, nonspin- 
ning BH binaries embedded in gas clouds under different assumptions for the 
motion of the binary with respect to the gas. They found evidence that 
the accretion rate and luminosity due to bremsstrahlung and synchrotron 
emission would be enhanced with respect to a sing le B H of the same mass 
as the binary, possibly being detectable by LSST |295| . Recent work sug- 



gest that the circumbinary disk surroundi ng BH bi naries should not produce 



detectable electromagnetic counterparts |294l . |296| . However, the magnetic 
field s produced by the circumbinary disk could affect the binary dynamics 
produce stronger electromagnetic counterparts 
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298 and extract en- 



ergy from the orbiting BHs, which ultimately merge within the standard 
Blandford-Znajek scenario, generating electromagnetic emission along dual 
or single jets that coul d be observable to large distance [299|, l300l . l30ll 1302 
(but see also (ioJ-IJO^l 

For further reading, we recommend the following reviews. An excellent 
summary of BH mass and spin measurements and (possible) evidence for 
event hor izon s based on "traditional" electromagnetic astronomy is given by 
Narayan 305| . The use of electromagnetic observations of BH s and neutron 



stars to probe strong-field gravity is reviewed by Psaltis |306| . A more ex- 
tended summary of nu merical results on gravitational recoil can be found in 
Zlochower et al. 285| . For a more g eneral rey iew of numerical BH simula- 



tions we recommend Centrella et al. 190, 307 



5. Black holes and high energy physics 

The breakthroughs in numerical relativity have opened the door to tackle 
many fundamental problems in BH physics and to address questions of wider 
interest. A new Golden Era in BH physics is just starting. We summarize 
below what we think have been the main developments in this relatively short 
period of time. 
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Mathematical physics and fundamental issues. A decades-old prob- 
lem in GR concerns the high-energy collision of two BHs. Due to the domi- 
nance of the gravitational interaction at large energies, this process is thought 
to describe general trans-planckian scattering of particles: at very large center 
of mass (CM) energies all interactions are "frozen" while gravity is boosted, 
thus all the details about internal structure of the colliding particles should be 
washed out. This is further supported by Thome's hoop conjecture (recently 



tested in highly dynamical situations |73|) , wh ich predicts BH formation from 



generic high-energy collisions of particles |308j . Because the final BH horizon 



cloaks all details about the structure of the colliding objects, matter does not 
matter and one can for simplicity study BHs as representing a wide class of 



colliding objects [38|, |39 



For large CM energies, BH collisions are arguably the most violent and 
nonlinear process one could conceive of. Evolving Einstein's equations in 
such extreme regimes poses new problems, such as the need to deal with all 
the different scales involved, and raises several questions. A fundamental 
issue concerns Cosmic Censorship. In four spacetime dimensions, BHs have 
an upper bound on their angular momentum, given bjo 

X^Sc/iGM')<l. (17) 

Is it possible to tune the impact parameter in such a way as to produce a final 
object spinning above the Kerr bound, or does nature somehow conspire to 
always radiate enough angular mom entum as to produc e a BH? This problem 
has been addressed in recent years 309|, ISlOl . Isill . |312 . 



We show the trajectory of two equal-mass, nonspinning, colliding BHs in 
Fig. m for CM velocities of around 0.75c. We notice three distinct regimes. 
For impact parameters above the scattering threshold {b > 6scat, the gravita- 
tional interaction between the two colliding BHs is relatively weak and they 
scatter to infinity. For impact parameters below the threshold of immediate 
merger {b < 6*), the BHs promptly collide and merge into a single BH. The 
transition between these two states for b^ < b < 6scat is a nonprompt merger, 
which puts the BHs in a so-called zoom-whirl orbit, i.e. the number of orbits 
rzorb ~ In 1 6—6* I . The results show that delayed mergers are necessary in order 



^We note that the Earth's spin XEarth ^ 10'^, while that of a spinning top Xtop ^ 10^^, 
so (even though the event horizon generators move at the speed of hght) BHs rotate 
"slowly" , as measured by their Kerr parameters. 
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Figure 2: Puncture trajectories of one BH for a scattering orbit (6[= 3.40M] > 6scat), a 
prompt merger (&[= 3.34Af] < b*) and a nonprompt merger (&* < b[= 3.39A'/] < fescat)- 
Taken from Ref. 131011. 
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Figure 3: Total energy radiated (left) and final BH spin (right) vs. impact parameter from 
sequence 2, the latter calculated using several methods. The vertical dashed green (dotted 
red) line is the estimated threshold of immediate merger b* (the scattering threshold &scat)- 



for the system to radiate excess angular momentum, so that the end product 
obeys flT7|) . The radiated energy and angular momentum increase strongly 
as the impact parameter approaches the threshold of immediate merger, as 
can be seen in Fig. [Hi For head-on collisions, numerical results suggest that 
a fracti on 1 4% of the CM energy is radiated in the limit of very large CM 



energy 



3091 



For finely tuned impact parameters, the total radiated energy 
(in GWs) can be of order 30% for a collision with v ~ 0.75c in the CM frame, 
while the spin of the final hole can reach ~ 96% of the bound (II 7p . An open 
issue is how much CM energy can be radiated for larger boosts. Finally, 
peak luminosities of around O.lc^/G are attained in these simulations. This 
number is orders of magnitude above any electro magnetic phenomena and 



close to the maximum conjectured possible value |313l . |314| |. The scattering 



threshold as a function of boost was studied by Shibata et al. and Sperhake 
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et al. |310l I3III |312| . For the range of boosts studied it is well approximated 

by 

fescat/M~(2.5±0.05)/t;. (18) 

Here M is the total ADM mass and v the velocity of each hole as measured 
in the CM frame. 

Higher-dimensional black holes. Although BH solutions and their prop- 
erties should in principle fall under the "mathematical physics" category, the 
activity in this specific field has been so intense that it deserves a subsec- 
tion of its own. Asymptotically fiat higher-dimensional black objects have 
a much richer structure than their four- dimensional counterparts. For in- 
stance, spherical topology is not the only allowed topology for objects with a 
horizon. One can also have, e.g., black rings, with a do nut-like topology. Re- 
markably, th ese t wo different horizon topologies coexist for certain regions in 



phase-space (315|. Explorations of the stability of general higher- dimensional 
BHs are in their infancy. Generically it has been co njec tured that, for D > 6, 
ultra-spinning Myers- Perry BHs will be unstable 316| . This instability has 
been confirrned by an analysis of linearized axisymmetric perturbations in 
D = 7,8,9 317|. Clearly, the study of the nonlinear development of these 
instabilities requires numerical methods, such as the ones reviewed here. A 
study of this type was very recently presented for nonaxisymmetric pertur- 



bations in D = 5 [318|, |319| , where it was found that a single spinning five 



dimensional Myers- Perry BH is unstable, for sufficientl y large rotation pa- 



rameter (confirming previous conjectures 320l . l32ll . |322 



General equilibrium states in anti-de Sitter backgrounds only recently 
started to be explored: nonaxisymmetric solutions have fin ally been built 



323|, confirming previous conjectures about their existence |32ll |322| . and 



Fully dynamical 



braneworld BHs are now also starting to be explored |324 . 
situations are still uncharted territory. 

Te V-scale gravity scenarios. The above studies are of direct relevance to 
some high-energy scenarios. An outstanding problem in high-energy physics 
is the extremely large ratio between the four dimensional Planck scale, 10^^ 
GeV, and the electroweak scale, 10^ GeV. It has been proposed that this 
hierarchy problem can be resolved by c onfining the Standard Model to a 



brane in a higher dimensional space (36|, l325l 1371 . |326 . 

In such models, the fundamental Planck scale - the energy at which 
gravitational interactions become strong - could be as low as 1 TeV. Thus, 
high-energy colliders, such as the Large Hadron Collider (LHC), may directly 
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probe strongly coupled gravitational physics |327l . l4ll . l38l |39| . In fact, such 
tests may even be routinely available in t he collisions of ultra high-energy 
cosmic rays with the Earth's a t mos phere 328l . l329l . |330| . or in astrophysi- 
cal BH environments 33ll l332l |333| (for reviews see J334l . l335l . |336|). The 
production of BHs at trans-Planckian collision energies (compared to the 
fundamental Planck scale) should be well described by using classical GR 



extended to D dimensions |41|, |38|, 139|, 1334 I335L l336j . The challenge is then 



to use the classical framework to determine the cross section for production 
and, for each initial setup, the fractions of the collision energy and angular 
momentum that are lost in the higher- dimensional space by GW emission. 
This information will be of paramount importance to improve the modelling 
of microscopic BH production in event generators. 

The first models for BH production in parton-parton collisions used a sim- 
ple black disk approach to estimate the cross section for production 38|, |39 . 
As we already described, only recently exact results for highly relativistic 
collisio ns where obtained i n four dimensions, using numerical relativity tech- 

The e xtension to hig her-dimensional spacetimes 

The first numerical results 



309, 310, 311, 312 



niques 

is a topic of current investigations |337l l338l . 1339 



339. 340 



ini- 



concerning low-energy collisions have been reported last year 

tial data for boosted E IH b i nari es have been constructed using an extended 

puncture approach in 34ll. |342|| and some results for high-energy collisions 



have been presented in [343 



AdS/CFT and holography. In 1997-98, a powerful new technique known 
as the AdS/CFT correspondence or, more generally, the gauge-gravity dual- 



ity, was introduced and rapidly developed [33[. This holographic correspon- 
dence provides an effective description of a nonperturbative, strongly coupled 
regime of certain gauge theories in terms of higher-dimensional classical grav- 
ity in AdS backgrounds. In particular, equilibrium and nonequilibrium prop- 
erties of strongly coupled thermal gauge theories are related to the physics of 
higher-dimensional BHs, black branes and their fluctuations. These studies 
revealed intriguing connections between the dynamics of BH horizons and 
hydrodynamics, and offer new perspectives on notoriously difficult problems, 
such as the BH information loss paradox, the nature of BH singularities and 
quantum gravity. 

Numerical relativity in anti-de Sitter backgrounds is bound to contribute 
enormously to our understanding of the gauge-gravity duality and is l ikely 



to have important applications in the interpretation of observations [35|, |344J . 
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345l |346| ■ For instance, in the context of the gauge-gravity duahty, high- 
energy colhsions of BHs have a dual description in terms of a) high-energy 
colhsions with balls of deconfined plasma surrounded by a confining phase, 
and b) the rapid localized heating of a deconfined plasma. These are the 
type of events that may have direct observational consequences f or the ex - 
periments at Brookhaven's Relativistic Heavy Ion Collider (RHIC) |345l . 1346 . 
Numerical relativity in anti-de Sitter is par ticularly difficult, and so far only 



very special situations have been handled 3471 . l348l . l349l l350l . l35ll |. Among 



these, we note the rece nt numeri cal proof that pure AdS space is unstable 
against collapse to BHs |353l . 13541 . 

BHs i n the context of high-energy physi cs have also been discussed by 



Pretorius 188 



356 . A recent review of numerical 



and in the white paper 
simulations of black strings in five spacetime dimensions is given by Lehner 
and Pretorius 



359 



dimensio ns is given in |357l . 1358 



the LHC 336 



An overview about numerical results on BHs in D > 4 
Finally we note Kanti's review on BHs at 



6. Conclusions 

We conclude this review with a brief summary of the most urgent prob- 
lems to be tackled by numerical relativity in the three main areas in the near 
future. 

Gravitational wave physics: In order to meet tight accuracy require- 
ments in GW template generation, especially for parameter estimation, we 
either need longer numerical waveforms or a denser spacing of templates in 
the parameter space. Both approaches will require computational resources 
that grow non-linearly as a function of the accuracy, either becaus e of the 
slow nature of the inspiral at larger BH separations (see Sec. 2 in 182| for 
a quantitative discussion) or because of the dimensionality of the parameter 
space. A second major task is the extension of existing template models to 
generic binaries with precessing spins and/or smaller mass ratios, bridging 
the gap to perturbative modeling of extreme mass ratio binaries. 

Astrophysics: While existing formulae for the kick and final spin resulting 
from the coalescence of BHs have been helpful for astroph ysical stu dies, the 

requires 



360, 361 



calibration of generic models as proposed by Boyle et al. 

a more comprehensive exploration of the parameter space. The relatively 

young field of numerical relativity simulations of BHs surrounded by accretion 
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disks will likely improve our insight into expected optical counterparts to BH 
binary mergers, and thus provide a vital tool for the area of multi-messanger 
astrophysics. 

High-energy physics: Arguably the most urgent task is the determination 
of the scattering cross-section and the loss of energy and momentum in GWs 
in trans-Planckian scattering of BHs in arbitrary dimensions. In view of the 
stability problems encountered in present studies, it also appears desirable 
to obtain a better understanding of the well-posedness of the formulations 
currently employed for such studies. Second, the modeling of BHs in generic 
spacetimes such as de Sitter and anti-de Sitter is still in its infancy, even 
if preliminary studies of BHs in non-asymptotically flat spacetimes are en- 
couraging. A better understanding of boundary issues, well-posedness of the 
formulations and diagnostics such as wave extraction tools and horizon find- 
ers will be required to open up this uncharted and fertile ground of research. 
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